ftp.cs.arizona.edu/ftp.cs.arizona.edu.tar/ftp.cs.arizona.edu/icon/newsgrp/group97b.txt/000049_icon-group-sender _Tue Aug 5 09:41:32 1997.msg

Internet Message Format | 2000-09-20 | 20KB

Received: from kingfisher.CS.Arizona.EDU by cheltenham.cs.arizona.edu; Tue, 5 Aug 1997 12:21:53 MST
Received: by kingfisher.CS.Arizona.EDU (5.65v4.0/1.1.8.2/08Nov94-0446PM)
    id AA15282; Tue, 5 Aug 1997 12:21:52 -0700
Message-Id: <[e-mail]>
X-Mailer: Novell GroupWise 4.1
Date: Tue, 05 Aug 1997 09:41:32 -0500
From: Charlie Hethcoat <[e-mail]>
To: [e-mail]
Cc: [e-mail]
Subject: Icon Engineering examples? -Reply
Mime-Version: 1.0
Content-Type: text/plain
Content-Disposition: inline
Errors-To: [e-mail]
Status: RO

I meant to answer this a couple of weeks ago.  Sorry about the delay.

I am a Senior Engineer in the Analysis group at Oceaneering Space
Systems.  We do a lot of stress and thermal analysis, and are heavily into
finite element modelling of astronaut crew equipment.  While my job isn't
primarily to program stuff, it is an activity I enjoy, and I occasionally need to
do some ad-hoc things with Icon.

As one example, here is a program that computes 2-dimensional area
properties of polygons of any size, with any number of embedded round
holes.  The program is fairly common (I first saw it programmed for a HP-67
calculator!), but nobody seems to think it should be a standard part of
anybody's library.  Hence, when I need it, I always have to reimplement it.

In my Icon version, there is a nice ascii syntax for the input, it handles any
number of vertices and holes, prints a neat little report, and generates a
gnuplot input file and control file for plotting the shape.  Of course, you'll
need gnuplot to plot it, or some plotting program that can handle gnuplot
input.  There are some weaknesses in the input code that I need to improve
(it isn't totally bulletproof), but the problem input syntax is so simple that
there is almost no risk of screwing it up.

Some people here complain that it isn't a CAD program, and they would
prefer to use, e. g. Autocad or Pro-Engineer.  Fine, let 'em, but I can be in
and out of my little old Icon program, with the answers,  in less time than it
takes them to log into Pro-E.  'Course, my program is only 2-D.  If 3-D is
needed, you have to go elsewhere.  But then the descriptive geometry
overhead becomes much more severe.

As you can see, Icon is a very nice language for general purpose mucking
about with numerics.  The only problem I have noticed is a tendency of the
left(), right(), and center() functions to truncate (rather than round properly)
floating point numbers when printing them.  Internal accuracy, though, seems
to be correct.

Warning!  This program requires command-line skills!!!   8^)  Maybe
someone can modify it to use vib--my skills to do that are yet insufficient.  I'd
like to see how my code can be adapted to that environment.
_____________________________________________

# aprop.icn - Compute and plot area properties of 2D polygons w/holes

link "w:/software/icon/mylib/scantool"     # Update paths for your system
link "w:/software/icon/mylib/numbtool"
link "w:/software/icon/mylib/plotfns"
link "basename"  # Std. Icon library
link "real2int"     # for floor() and ceil()
link "w:/software/icon/mylib/minmax"    # for minlist() and maxlist()

record holetype(x, y, d) # each hole gets location (x, y) and diameter d

global infile, outfile        # Text input and output
global plotname, plotfile     # Plot data output
global gnuname, gnufile       # Plot control output for Gnuplot
global bname  # For access by plotting routines
global x, y, holelist  # lists
global i, n
global xmin, xmax, xrange, ymin, ymax, yrange # For plotting ranges

global area, xcg, ycg
global Ixx, Iyy, Ixy
global Ixcg, Iycg, Ixycg
global angle
global Ixprime, Iyprime, Ixyprime
global J

procedure main(args)

     static rcsid
     initial {
          rcsid := "$Id: aprop.icn 1.4 1996/12/18 21:16:39 Hethcoat Exp
Hethcoat $"
     }

     #
     # $Log: aprop.icn $
     # Revision 1.4  1996/12/18 21:16:39  Hethcoat
     # Noticed a typo in the name of the rcsid variable.
     #
     # Revision 1.3  1996/10/10 22:35:27  Hethcoat
     # Plotting enhancements:  Now generates all gnuplot data,
     # scale factors for plotting now fully data-driven and are always in
     # multiples of 0.5 world units per inch.
     #
     # Revision 1.2  1996/10/05 21:35:29  Hethcoat
     # Merely added the RCS keyword strings.
     #
     #

     infile := &input
     outfile := &output

     if *args = 0 then {  # A kluge for you to enjoy...
          system("list.com aprop.hlp")
          exit(0)
     }

     # Looks for an input file in the form modelname.ext
     # Generates output files in the same format.
     bname := basename(args[1], ".")
     infile := open(args[1], "r") | stop("Cannot open ", args[1], ".")
     outname := bname || ".apo"
     plotname := bname || ".dat"
     gnuname := bname || ".gnu"
     outfile := open(outname, "w") | stop("Cannot open ", basename, ".")

     # Setup

     holelist := []  # add holes as found
     x := []
     y := []
     i := 0
     write(outfile, center("Input Data Listing", 72, "-"))
     write(outfile, )
     # Input phase
     while line := decomm(infile, outfile) do {
          line ? {
               if tab(match("hole")) then {
                    ws(); xhole := num()
                    ws(); yhole := num()
                    ws(); dhole := num()
                    put(holelist, holetype(xhole, yhole, dhole))
                    next
               }

               i +:= 1
               put(x, num())
               ws()
               put(y, num())
          }
     }

     n := *x  # The number of points read, including any repeats.

     # Need ranges of the input data for plotting purposes:
     xmin   := minlist(x)
     xmax   := maxlist(x)
     xrange := xmax - xmin
     ymin   := minlist(y)
     ymax   := maxlist(y)
     yrange := ymax - ymin

     # Computations
     do_numbers()

     # Report
     do_report()

     # Plot control file
     outplot1(xmin, xmax, ymin, ymax)

     # Data plot output
     outplot2()

     exit(0)

end

# do_numbers - compute and print the detailed stuff
procedure do_numbers()
     area := xcg := ycg := Ixx := Iyy := Ixy := 0.0
     n := *x
     every i := (1 to n-1) do {
          delta_x := x[i+1] - x[i]
          sum_x   := x[i+1] + x[i]
          delta_y := y[i+1] - y[i]
          sum_y   := y[i+1] + y[i]
          area +:= delta_y*sum_x
          xcg +:= delta_y*(sum_x^2 + delta_x^2/3.)
          ycg +:= delta_x*(sum_y^2 + delta_y^2/3.)
          Ixx +:= delta_x*sum_y*(sum_y^2 + delta_y^2)
          Iyy +:= delta_y*sum_x*(sum_x^2 + delta_x^2)
          # Ixy is so big it has to be done in pieces:
          term1 := x[i+1]*y[i] - x[i]*y[i+1]
          term2 := x[i+1]^2 + x[i]^2
          Ixy1 := delta_y^2*sum_x*(term2)/8.
          Ixy2 := delta_y*(term1)*(term2 + x[i+1]*x[i] )/3.
          Ixy3 := (term1)^2*sum_x/4.
          if not wee(delta_x) then
               Ixy +:= (Ixy1 + Ixy2 + Ixy3)/delta_x
          # else there is no contribution when delta_x
          # is zero.
     }
     area /:= -2.0
     xcg /:= (-8.*area)
     ycg /:= (8.*area)
     Ixx /:= 24.
     Iyy /:= (-24.)
     # Ixy is OK as it comes out.

     # The above is correct before adding round holes.
     # Now have to correct for holes found, if any:

     every hole := !holelist do {
          holearea := -&pi*hole.d^2/4.0 # Note:  negative area for holes!!
          holeIxx := -&pi*hole.d^4/64.  # also negative!!
          holeIyy := holeIxx # by symmetry
          holeIxy := 0.0  # also by symmetry
          newarea := area + holearea # don't tamper with area just yet!
          xcg := (xcg*area + holearea*hole.x)/newarea # updated x c.g.
          ycg := (ycg*area + holearea*hole.y)/newarea
          Ixx +:= (holeIxx + holearea*hole.y^2)
          Iyy +:= (holeIyy + holearea*hole.x^2)
          Ixy +:= (holeIxy + holearea*hole.x*hole.y)
          area := newarea
     }

     # Finish up with the final numbers:

     Ixcg := Ixx - area*ycg^2
     Iycg := Iyy - area*xcg^2
     Ixycg := Ixy - area*xcg*ycg
     angle := 0.5*atan(-2.*Ixycg, Ixcg - Iycg)
     c := cos(angle)
     s := sin(angle)
     s2 := sin(2.*angle)
     c2 := cos(2.*angle)
     Ixprime := Ixcg*c^2 + Iycg*s^2 - Ixycg*s2
     Iyprime := Iycg*c^2 + Ixcg*s^2 + Ixycg*s2
     J := Ixprime + Iyprime
     Ixyprime := 0.5*(Ixcg - Iycg)*s2 + Ixycg*c2
     return
end

# do_report - print the report
procedure do_report()
     write(outfile, )
     write(outfile, center("Area Properties Report", 72, "-"))
     write(outfile, )
     write(outfile, "Polygon:  ", n - 1, " points")
     write(outfile, "Holes:  ", *holelist)
     write(outfile, "area = ", area)
     write(outfile, "xcg  = ", xcg)
     write(outfile, "ycg  = ", ycg)
     write(outfile, "Original axes:")
     write(outfile, "Ixx  = ", Ixx)
     write(outfile, "Iyy  = ", Iyy)
     write(outfile, "Ixy  = ", Ixy)
     write(outfile, "Centroidal axes (parallel to Original):")
     write(outfile, "Ixxcg = ", Ixcg)
     write(outfile, "Iyycg = ", Iycg)
     write(outfile, "Ixycg = ",Ixycg)
     write(outfile, "Principal axes (rotated from Centroidal):")
     write(outfile, "angle = ", rtod(angle), " degrees")
     write(outfile, "Ix'x' = ", Ixprime)
     write(outfile, "Iy'y' = ", Iyprime)
     write(outfile, "Ix'y' = ", Ixyprime)
     write(outfile, "Polar Moment of Inertia:")
     write(outfile, "J = ", J)
     return
end

# outplot1 - Output file (*.gnu) with plotting controls
procedure outplot1(xmin, xmax, ymin, ymax)

     # Right now, this plots only to an hpljii at 150x150.  We'll
     # improve it later.

     static wide, high, clear

     local xrange, yrange
     local xbox, ybox
     local xmini, xmaxi, ymini, ymaxi
     local xsizei, ysizei
     local xscale, yscale
     local scale
     local max_x, max_y

     initial {
          # Need the actual plot area.
          # These are correct for the HP Laserjet 4 Plus
          # in portrait mode (so far!)
          wide := 5.35
          high := 6.06
          clear := 0.5  # minimum clearance around the data
     }

     # Imagine a box just enclosing the data:
     xbox := xmax - xmin
     ybox := ymax - ymin

     # Now find the next larger box with at least the stipulated clearance
     # and with integer dimensions:
     xmini := ceil(xmin - clear)
     xmaxi := ceil(xmax + clear)
     ymini := ceil(ymin - clear)
     ymaxi := ceil(ymax + clear)

     # Get the size of that integer box:
     xsizei := xmaxi - xmini
     ysizei := ymaxi - ymini

     # Compute separate x and y scale factors for plotting:
     xscale := xsizei/wide
     yscale := ysizei/high

     # The plot will be made to equal x and y scales, so it has to be
     # the larger of the two scales.  Furthermore, it would be convenient
     # for the scale to be a nice number, i. e. positive multiples of 0.5.
     scale := 0.5*ceil(2.0*maxlist([xscale, yscale]))

     # Now, get the plotting parameters in shape for use with the gnuplot
     # xrange and yrange commands:

     max_x := xmini + scale*wide
     max_y := ymini + scale*high

     gnufile := open(gnuname, "w") | stop("Cannot open ", gnuname, ".")
     write(gnufile,"set data style lines")
     write(gnufile,"set xrange [", xmini, ":", max_x, "]")
     write(gnufile,"set yrange [", ymini, ":", max_y, "]")
     write(gnufile,"set terminal hpljii 150")
     write(gnufile,"set output '", bname || ".prn'")
     write(gnufile,"plot '", plotname, " '")
     write(gnufile,"exit")
#    write(gnufile,"") PATTERN--KEEP!
     close(gnufile)
     return
end

# outplot2 - Output file (*.dat) of results
procedure outplot2()
     local length

     plotfile := open(plotname, "w") | stop("Cannot open " || plotname || ".")
     length := 1.25*minlist([xrange, yrange])

     # Plot the input (x, y) data:
     comment("Plot file generated by aprop on " || &dateline)
     comment("Cross section:  " || bname)
     comment("")  # empty comment
     comment("Polygon data:")
     drawdata(x, y)
     comment("Centroidal axes:")
     drawcross(xcg, ycg, length, 0.0)
     comment("Principal axes:")
     drawcross(xcg, ycg, length, angle)
     comment("End of data.")
     close(plotfile)
     return
end

# drawdata - draw (x, y) list data
procedure drawdata(x, y)
     local i, hole
     move(x[1], y[1])
     every i := 2 to *x do
          draw(x[i], y[i])
     every hole := !holelist do {
          comment("Hole:")
          drawcircle(hole.x, hole.y, 0.5*hole.d)
          comment("Crosshairs for preceding hole:")
          drawcross(hole.x, hole.y, hole.d, 0.0)
     }
     return
end
________________________________________

# numbtool.icn - basic numerical math aids (nothing too fancy)
# Usage:  line w:\software\icon\mylib\numbtool
# Author:      Charles L. Hethcoat III
# Date:        September 14, 1996

# macheps - discern machine epsilon by test
procedure macheps()
     static eps
     initial eps := 1.0
     while 1.0 + eps > 1.0 do
          eps *:= 0.5
     # Comes out only when eps is indiscernible, i.e. makes
     # no difference when added to 1.0.  This is more or less the
     # definition of the machine epsilon.
     return eps
end

# wee - succeed if number is too "wee" (love that word) to divide by
procedure wee(x)

     # Typical usage:  if not wee(x) then y := 1.0/x
     # (wee() has to succeed, so that not wee fails, if x is wee.)
     static thresh
     initial thresh := 5.0*macheps()
     if abs(x) < thresh then
          return
     # else fail
end
______________________________________

# scantool.icn - simple token-scanning tools for use in various projects
# Usage:  link w:\software\icon\mylib\scantool
# Author:      Charles L. Hethcoat III
# Date:        September 14, 1996

# id - scan an identifier
procedure id()
     static first, rest
     initial {
          first := &letters ++ '_'      # First chars in an identifier
          rest := first ++ &digits      # Remaining ones
     }
     return scantok(first, rest)
end

# num - scan and return the next number in &subject
procedure num()
     static first, rest
     initial {
          first := &digits ++ '+-.'    # First chars in a number
          rest := first ++ 'e'        # Remaining ones
     }
     return real(scantok(first, rest))
end

# scantok - scan tokens starting with first, finishing with rest
procedure scantok(first, rest)
     ws()
     return tab(any(first)) || tab(many(rest)) | tab(any(first))
end

# ws - skip white space (blanks and tabs)
procedure ws()
     tab(many(' \t'))
end

# decomm - echo lines from infile and return decommented ones only
procedure decomm(infile, outfile)
     local line
     while line := write(outfile, read(infile)) do
          line ? {
               ws()
               if ='#' | pos(0) then next
               return trim(tab(upto('#') | 0))
          }
end
_____________________________________________________
# minmax.icn:
############################################################################
#
#     These procedures return the maximum and minimum of a list containing
#  any number of numeric arguments.  Modified by CLH from the library
#  routine minmax, which uses variable-length argument lists.
#
############################################################################

procedure minlist(L)
   local minimum, n
   minimum := L[1] | fail
   every n := 2 to *L do
      minimum >:= L[n]
   return minimum
end

procedure maxlist(L)
   local maximum, n
   maximum := L[1] | fail
   every n := 2 to *L do
      maximum <:= L[n]
   return maximum
end
_____________________________________
# plotfns.icn:
# Plot routines that I can use in future programs.
# All of the following "plot routines" just generate Gnuplot output!!

# drawcross - plot a cross with bar length L centered at (x, y) tilted theta
procedure drawcross(x, y, L, theta)
     centerline(x, y, L, theta)              # line
     centerline(x, y, L, theta + 0.5*&pi)    # same line rotated +90 degrees
     return
end

# centerline - plot a line of length L centered at (x, y), tilted theta
procedure centerline(x, y, L, theta)
     local L2, dx, dy
     L2 := 0.5*L
     dx := L2*cos(theta)
     dy := L2*sin(theta)
     move(x - dx, y - dy)
     draw(x, y)
     draw(x + dx, y + dy)
     return
end

# drawcircle - plot a circle of radius r centered at (x, y)
procedure drawcircle(xc, yc, r)
     static n, deltheta
     local i, theta, x, y, x1, y1
     initial {
          n := 16    # Number of points to draw
          deltheta := 2.0*&pi/n
     }
     every i := 1 to n do {
          theta := real(i - 1)*deltheta
          x := xc + r*cos(theta)
          y := yc + r*sin(theta)
          if i = 1 then {
               x1 := x; y1 := y  # Save 1st point for doing the last line.
               move(x, y)  # and don't drag the pen.
          }
          else
               draw(x, y)
     }
     draw(x1, y1) # Close the circle.
     return
end

# move - move to a new starting coordinate (x, y) for the next line
procedure move(x, y)
     write(plotfile) # Gnuplot starts a new line with a blank input line.
     write(plotfile, x, " ", y)
     return
end

# draw - draw a line to final coordinate (x, y)
procedure draw(x, y)
     write(plotfile, x, " ", y)
     return
end

# comment - add a gnuplot comment line to the data
procedure comment(line)
     write(plotfile, "# ", line)
     return
end

_____________________________________

# aprop.inp - a right triangle with a rt. triangular hole & 2 round holes

# Sample input data
# Note the comment convention:  From any '#' to the end of a line is ignored.
# You may use comments in all your input files for documentation purposes.

# Start the outer triangle.  Move clockwise for areas:

# (X,Y) coordinate pairs:

3.   1.
3.   7.
14.  7.
3.   1. # Closure on outer triangle--necessary!

# Polygonal holes are defined by moving to and fro along a "cut" line.

4.   4. # Move along cut to inner triangle

# Polygonal holes are done the same way, except moving
counterclockwise.

9.   6. # Starts the inner triangle
4.   6.
4.   4. # Closure on inner triangle
3.   1. # Move outward on cut to final closure

# Finally, let's have two round holes:

#    X         Y         Diameter

hole 5.5       3.5       1.0
hole 10.0      6.0       0.5
_______________________________

That's it.  Good luck.

Charles Hethcoat

>>> "Brian A. Rogers" <[e-mail]> 1997 Jul 24,
02:05pm >>>
 Hi...

I'm looking for some examples of Icon programs that prvide some sort of
engineering or science function.  From what I've read so far, it's a neat
language.  I'd like to see someone using it, though.


--
----------------------------------------------------------
Brian A. Rogers
Student Manager, Laptop Orientation
Rose-Hulman Institute of Technology
[phone]
-----------------------------------------------------------
"There are no problems, only solutions."
Tron, 1982
-----------------------------------------------------------